library(tidyverse)
library(lmtest)
library(sandwich)
library(stargazer) #for tables
library(marginaleffects)
library(gridExtra)
#Remember to set WD


##########################################################################################
#Government Satisfaction DV
#load data
df1 <- read.csv(("data\\clean\\roe2023data.csv"))

#set FR as base category
df1$prtvtfam <- relevel(factor(df1$prtvtfam), ref = 4) 

#no fixed effects, no controls
reg.1.basic <- lm(stfgov ~ coalition.current + prtvtfam + coalition.current*prtvtfam, 
                  data = df1)

summary(reg.1.basic)

#robust standard errors
se.reg.1.basic <- coeftest(reg.1.basic, vcov = vcovHC)
se.reg.1.basic


#controls, no fixed effects
reg.1.cont <- lm(stfgov ~ coalition.current + prtvtfam + coalition.current*prtvtfam + agea + blgetmg + eduyrs +
                   gndr + hinctnta + lrscale + mbtru + seat_share, 
                 data = df1)

summary(reg.1.cont)

#robust standard errors
se.reg.1.cont <- coeftest(reg.1.cont, vcov = vcovHC)
se.reg.1.cont




#country and year fixed effects, controls (full model)
#set FR as base category

reg.1.fixed <- lm(stfgov ~ coalition.current + prtvtfam + coalition.current*prtvtfam + agea + blgetmg + eduyrs +
                    gndr + hinctnta + lrscale + mbtru + seat_share + 
                    factor(cntry) + factor(inwyys), 
                  data = df1)

p <- predictions(reg.1.fixed, 
                 variables = list(coalition.current = 0:1),
                 newdata = datagrid(prtvtfam = unique))

#Figure 1 (first panel)
stfgovplot1 <- ggplot(p,
       aes(x = estimate,
           y = factor(prtvtfam,
                      level = c("Communist",  "Social Dem", 
                                "Liberal", "Christian Dem", 
                                "Conservative", "Far-Right")),
           xmin = conf.low,
           xmax = conf.high,
           color = factor(coalition.current))) +
  geom_errorbarh(height = 0,
                 position = position_dodge(width = .2),
                 show.legend = F,
                 size = 1) +
  geom_point(position = position_dodge(width = .2),
             show.legend = F,
             size = 2) +
  labs(y = "Party Family",
       x = "Predicted Government Satisfaction with 95% CI") +
  scale_color_grey(start = .6,
                   end = 0) + 
  scale_x_continuous(breaks = round(seq(4, 7, by = 0.5),1)) +
  theme_bw() 

c <- comparisons(reg.1.fixed, 
                 variables = list(prtvtfam = "reference"), #what to compare to each other, creates contrasts across values
                 newdata = datagrid(coalition.current = 0:1),
                 comparison = 	\(hi, lo) lo - hi,#what do you want contrasts for?
                 vcov = "HC3") 

#Figure 2 (first panel)
stfgovplot2 <- ggplot(c,
       aes(x = estimate,
           y = factor(contrast,
                      level = c("Communist, Far-Right",  "Social Dem, Far-Right", 
                                "Liberal, Far-Right", "Christian Dem, Far-Right", 
                                "Conservative, Far-Right")),
           xmin = conf.low,
           xmax = conf.high,
           color = factor(coalition.current))) +
  geom_errorbarh(height = 0,
                 position = position_dodge(width = .2),
                 show.legend = F,
                 size = 1) +
  geom_point(position = position_dodge(width = .2),
             show.legend = F,
             size = 2) +#change in satsifaction by moving from 1 party to another
  labs(y = "Party Family Contrasts",
       x = "Differences in Expected Government Satisfaction with 95% CI") +
  scale_color_grey(start = .6,
                   end = 0) + 
  theme_bw()


c <- comparisons(reg.1.fixed, 
                 variables = list(coalition.current = 0:1),
                 newdata = datagrid(prtvtfam = unique),
                 vcov = "HC3") %>%
  glimpse()

#Figure 3 (first panel)
stfgovplot3 <- ggplot(c,
       aes(x = estimate,
           y = factor(prtvtfam,
                      level = c("Communist",  "Social Dem", 
                                "Liberal", "Christian Dem", 
                                "Conservative", "Far-Right")),
           xmin = conf.low,
           xmax = conf.high)) +
  geom_errorbarh(height = 0,
                 size =1) +
  geom_point(size = 2) +
  labs(y = "Party Family",
       x = "Within-Party Difference in Expected Government Satisfaction with 95% CI") +
  theme_bw()



c <- comparisons(reg.1.fixed, 
                 variables = list(coalition.current = 0:1),
                 newdata = datagrid(prtvtfam = unique),
                 hypothesis = "revreference",
                 vcov = "HC3") %>%
  glimpse()

#Figure 4 (first panel)
stfgovplot4 <- ggplot(c,
       aes(x = estimate,
           y = factor(term,
                      level = c("Row 1 - Row 4",  "Row 1 - Row 6",
                                "Row 1 - Row 2", "Row 1 - Row 3",
                                "Row 1 - Row 5")),
           xmin = conf.low,
           xmax = conf.high)) +
  geom_errorbarh(height = 0,
                 size = 1) +
  geom_point(size = 2) +#how much more serving in govt boosts satisfaction
  labs(y = "Party Family Contrasts",
       x = "Difference in Expected Government Satisfaction with 95% CI") +
  scale_y_discrete(labels=c("Row 1 - Row 5" = "Far-right, Conservative", "Row 1 - Row 3" = "Far-Right, Christian Dem",
                            "Row 1 - Row 2" = "Far-Right, Liberal", "Row 1 - Row 6" = "Far-Right, Social Dem",
                            "Row 1 - Row 4" = "Far-Right, Communist")) +
  theme_bw()


summary(reg.1.fixed)

#robust standard errors
se.reg.1.fixed <- coeftest(reg.1.fixed, vcov = vcovHC)
se.reg.1.fixed


stargazer(se.reg.1.basic, se.reg.1.cont, se.reg.1.fixed,
          keep = c("coalition.current",
                   "prtvtfam",
                   "seat_share",
                   "Constant"),
          dep.var.labels = c("Satisfaction with Government"),
          covariate.labels = c("Government", "Christian Dem.", "Communist",
                               "Conservative", "Liberal", "Social Dem.", "Seat Share",
                               "Government $\times$ Christian Dem.", "Government $\times$ Communist ",
                               "Government $\times$ Conservative", "Government $\times$ Liberal", 
                               "Government $\times$ Social Dem.",
                               "Constant"),
          notes.label = "Robust standard errors in parentheses",
          add.lines = c("N", "Adj. R-squared", "F-statistic", "Control Variables","Fixed Effects"),
          table.placement = "H",
          font.size = "tiny")






##########################################################################################################
#Party Closeness DV

#load data
df1 <- read.csv(("data\\clean\\roe2023data.csv"))


#set FR as base category
df1$prtvtfam <- relevel(factor(df1$prtvtfam), ref = 4) 


#DV = probability that the party you feel close to is the same one you voted for
reg.2.basic <- glm(clsprty.same ~ f.coalition.current + factor(prtvtfam) + f.coalition.current*factor(prtvtfam),
                   data = df1,
                   family = binomial(link = "logit"))
summary(reg.2.basic)

#robust standard errors
se.reg.2.basic <- coeftest(reg.2.basic, vcov = vcovHC)
se.reg.2.basic
logLik(reg.2.basic)

reg.2.cont <- glm(clsprty.same ~ f.coalition.current + factor(prtvtfam) + f.coalition.current*factor(prtvtfam) +
                    agea + blgetmg + eduyrs +
                    gndr + hinctnta + lrscale + 
                    mbtru + seat_share,
                  data = df1,
                  family = binomial(link = "logit"))

summary(reg.2.cont)

#robust standard errors
se.reg.2.cont <- coeftest(reg.2.cont, vcov = vcovHC)
se.reg.2.cont
logLik(reg.2.cont)







reg.2.full <- glm(clsprty.same ~ f.coalition.current + factor(prtvtfam) + f.coalition.current*factor(prtvtfam) +
                    agea + blgetmg + eduyrs +
                    gndr + hinctnta + lrscale + 
                    mbtru + seat_share +
                    f.cntry + f.inwyys, 
                  data = df1,
                  family = binomial(link = "logit"))
summary(reg.2.full)
#robust standard errors
se.reg.2.full <- coeftest(reg.2.full, vcov = vcovHC)
se.reg.2.full
logLik(reg.2.full)

p <- predictions(reg.2.full, 
                 variables = list(f.coalition.current = 0:1),
                 newdata = datagrid(prtvtfam = unique)) %>%
  glimpse()

#Figure 1 (second panel)
closeplot1 <- ggplot(p,
       aes(x = estimate,
           y = factor(prtvtfam,
                      level = c("Communist",  "Social Dem", 
                                "Liberal", "Christian Dem", 
                                "Conservative", "Far-Right")),
           xmin = conf.low,
           xmax = conf.high,
           color = factor(f.coalition.current))) +
  geom_errorbarh(height = 0,
                 position = position_dodge(width = .2),
                 size = 1) +
  geom_point(position = position_dodge(width = .2),
             size = 2) +#need to dodge
  labs(y = "",
       x = "Predicted Party Closeness with 95% CI",
       color = "Executive Participation") +
  scale_color_grey(start = .6,
                   end = 0,
                   labels=c("Opposition", "In Government")) + 
  theme_bw()

#Figure 2 (second panel)
c <- comparisons(reg.2.full, 
                 variables = list(prtvtfam = "reference"), #what to compare to each other, creates contrasts across values
                 newdata = datagrid(f.coalition.current = 0:1),
                 comparison = 	\(hi, lo) lo - hi) %>% #what do you want contrasts for?
  glimpse()

closeplot2 <- ggplot(c,
       aes(x = estimate,
           y = factor(contrast,
                      level = c("Communist, Far-Right",  "Social Dem, Far-Right", 
                                "Liberal, Far-Right", "Christian Dem, Far-Right", 
                                "Conservative, Far-Right")),
           xmin = conf.low,
           xmax = conf.high,
           color = factor(f.coalition.current))) +
  geom_errorbarh(height = 0,
                 position = position_dodge(width = .2),
                 size = 1) +
  geom_point(position = position_dodge(width = .2),
             size = 2)+ #change in closeness by moving from 1 party to another
  labs(y = "",
       x = "Differences in Expected Party Closeness with 95% CI",
       color = "Executive Participation") +
  scale_color_grey(start = .6,
                   end = 0,
                   labels=c("Opposition", "In Government")) + 
  theme_bw()

#compares closeness among parties not serving in government


c <- comparisons(reg.2.full, 
                 variables = list(f.coalition.current = 0:1),
                 newdata = datagrid(prtvtfam = unique)) %>%
  glimpse()

#Figure 3 (second panel)
closeplot3 <- ggplot(c,
       aes(x = estimate,
           y = factor(prtvtfam,
                      level = c("Communist",  "Social Dem", 
                                "Liberal", "Christian Dem", 
                                "Conservative", "Far-Right")),
           xmin = conf.low,
           xmax = conf.high)) +
  geom_errorbarh(height = 0,
                 size = 1) +
  geom_point(size = 2) +
  labs(y = "",
       x = "Within-Party Difference in Expected Party Closeness with 95% CI") +
  theme_bw()



c <- comparisons(reg.2.full, 
                 variables = list(f.coalition.current = 0:1),
                 newdata = datagrid(prtvtfam = unique),
                 hypothesis = "revreference") %>%
  glimpse()

#figure 4 (second panel)
closeplot4 <- ggplot(c,
       aes(x = estimate,
           y = factor(term,
                      level = c(  "Row 1 - Row 3", "Row 1 - Row 6",
                                  "Row 1 - Row 5",
                                  "Row 1 - Row 2", "Row 1 - Row 4")),
           xmin = conf.low,
           xmax = conf.high)) +
  geom_errorbarh(height = 0,
                 size = 1) +
  geom_point(size = 2)+ #how much more serving in govt boosts satisfaction
  labs(y = "",
       x = "Difference in Expected Party Closeness with 95% CI") +
  scale_y_discrete(labels=c("Row 1 - Row 5" = "Far-right, Liberal", "Row 1 - Row 3" = "Far-Right, Communist",
                            "Row 1 - Row 2" = "Far-Right, Christian Democrat", "Row 1 - Row 6" = "Far-Right, Social Dem",
                            "Row 1 - Row 4" = "Far-Right, Conservative")) +
  theme_bw()




stargazer(se.reg.2.basic, se.reg.2.cont, se.reg.2.full,
          keep = c("f.coalition.current",
                   "prtvtfam",
                   "seat_share",
                   "govt_eff",
                   "Constant"),
          dep.var.labels = c("Party Closeness"),
          covariate.labels = c("Government", "Christian Dem.", "Communist", "Conservative",
                               "Liberal", "Social Dem.", "Seat Share", 
                               "Government $\times$ Christian Dem.", "Government $\times$ Communist", 
                               "Government $\times$ Conservative", "Government $\times$ Liberal", "Government $\times$ Social Dem.", 
                               "Constant"),
          notes.label = "Standard errors in parentheses",
          table.placement = "H")

pdf("figure1.pdf", width = 11.5, height = 4)
grid.arrange(stfgovplot1, closeplot1, ncol=2)
dev.off()

pdf("figure2.pdf", width = 12, height = 4)
grid.arrange(stfgovplot2, closeplot2, ncol=2)
dev.off()

pdf("figure3.pdf", width = 12.5, height = 4)
grid.arrange(stfgovplot3, closeplot3, ncol=2)
dev.off()

pdf("figure4.pdf", width = 11.5, height = 4)
grid.arrange(stfgovplot4, closeplot4, ncol=2)
dev.off()



#load data
df1 <- read.csv(("data\\clean\\roe2023data.csv"))


#set FR as base category
df1$prtvtfam <- relevel(factor(df1$prtvtfam), ref = 4) 
df1$f.prtvtfam <- relevel(factor(df1$f.prtvtfam), ref = 4) 

#ROBUSTNESS
#############################################################
#Triple interaction term
#country and year fixed effects, controls (full model)
reg.3.fixed <- lm(stfgov ~ factor(coalition.current) + factor(prtvtfam) + months.last.elec+ factor(coalition.current)*factor(prtvtfam)*months.last.elec + 
                    agea + blgetmg + eduyrs +
                    gndr + hinctnta + lrscale + mbtru + seat_share + 
                    factor(cntry) + factor(inwyys), 
                  data = df1)
summary(reg.3.fixed)

#robust standard errors
se.reg.3.fixed <- coeftest(reg.3.fixed, vcov = vcovHC)
se.reg.3.fixed

#probability that the party you feel close to is the same one you voted for
reg4.full <- glm(clsprty.same ~ f.coalition.current + f.prtvtfam + months.last.elec +
                   f.coalition.current*f.prtvtfam*months.last.elec +
                   agea + blgetmg + eduyrs +
                   gndr + hinctnta + lrscale + 
                   mbtru + seat_share +
                   f.cntry + f.inwyys, 
                 data = df1,
                 family = binomial(link = "logit"))
summary(reg4.full)

stargazer(se.reg.3.fixed, reg4.full,
          keep = c("coalition.current",
                   "prtvtfam",
                   "months.last.elec",
                   "Constant"),
          dep.var.labels = c("Satisfaction with Government", "Party Closeness"),
#          covariate.labels = c("Coalition", "Christian Dem.", "Liberal",
#                               "Social Dem.", "Communist", "Conservative", "Coalition1", "Christian Dem.1", "Liberal1",
#                               "Social Dem.1", "Communist1", "Conservative1", "Months", "Seat Share", "Govt Effectiveness",
#                               "Coalition $\times$ Christian Dem.", "Coalition $\times$ Liberal",
#                               "Coalition $\times$ Social Dem.", "Coalition $\times$ Communist", 
#                               "Coalition $\times$ Conservative", "Coalition $\times$ Months", "Coalition $\times$ Months",
#                               "Christian Dem $\times$ Months", "Liberal $\times$ Months", "Social Dem. $\times$ Months",
#                               "Communist $\times$ Months", "Conservative $\times$ Months"),
          notes.label = "",
          add.lines = c("N", "Adj. R-squared", "F-statistic", 
                        "Log-Liklihood", "Control Variables","Fixed Effects"),
          table.placement = "H",
          font.size = "tiny")



##############################################################
#exclluding specific countries
#NO SWITZERLAND
df.no.ch <- df1 %>%
  filter(cntry != "CH")

#country and year fixed effects, controls (full model)
reg.no.ch <- lm(stfgov ~ factor(coalition.current) + factor(prtvtfam) + factor(coalition.current)*factor(prtvtfam) + agea + blgetmg + eduyrs +
                  gndr + hinctnta + lrscale + mbtru  + seat_share + 
                  factor(cntry) + factor(inwyys), 
                data = df.no.ch)
summary(reg.no.ch)

#robust standard errors
se.reg.no.ch <- coeftest(reg.no.ch, vcov = vcovHC)
se.reg.no.ch

#NO HUNGARY OR POLAND
df.no.hupl <- df1 %>%
  filter(cntry != "HU") %>%
  filter(cntry != "PL")

#country and year fixed effects, controls (full model)
reg.no.hupl <- lm(stfgov ~ factor(coalition.current) + factor(prtvtfam) + factor(coalition.current)*factor(prtvtfam) + agea + blgetmg + eduyrs +
                    gndr + hinctnta + lrscale + mbtru + seat_share  + 
                    factor(cntry) + factor(inwyys), 
                  data = df.no.hupl)
summary(reg.no.hupl)

#robust standard errors
se.reg.no.hupl <- coeftest(reg.no.hupl, vcov = vcovHC)
se.reg.no.hupl

#NO SWITZERLAND, HUNGARY, OR POLAND
df.no.all <- df1 %>%
  filter(cntry != "HU") %>%
  filter(cntry != "PL") %>%
  filter(cntry != "CH")

#country and year fixed effects, controls (full model)
reg.no.all <- lm(stfgov ~ factor(coalition.current) + factor(prtvtfam) + factor(coalition.current)*factor(prtvtfam) + agea + blgetmg + eduyrs +
                   gndr + hinctnta + lrscale + mbtru + seat_share + 
                   factor(cntry) + factor(inwyys), 
                 data = df.no.all)
summary(reg.no.all)



#robust standard errors
se.reg.no.all <- coeftest(reg.no.all, vcov = vcovHC)
se.reg.no.all

stargazer(se.reg.no.ch, se.reg.no.hupl, se.reg.no.all,
          keep = c("coalition.current",
                   "prtvtfam",
                   "seat_share",
                   "govt_eff",
                   "Constant"),
          dep.var.labels = c("Satisfaction with Government"),
#          covariate.labels = c("Coalition", "Christian Dem.", "Liberal",
#                               "Social Dem.", "Communist", "Conservative", "Seat Share", "Govt Effectiveness",
#                               "Coalition $\times$ Christian Dem.", "Coalition $\times$ Liberal",
#                               "Coalition $\times$ Social Dem.", "Coalition $\times$ Communist", 
#                               "Coalition $\times$ Conservative",
#                               "Constant"),
          notes.label = "Robust standard errors in parentheses",
          add.lines = c("N", "Adj. R-squared", "F-statistic", "Control Variables","Fixed Effects"),
          table.placement = "H",
          font.size = "tiny")




#probability that the party you feel close to is the same one you voted for
reg2.no.ch <- glm(clsprty.same ~ f.coalition.current + f.prtvtfam + f.coalition.current*f.prtvtfam +
                    agea + blgetmg + eduyrs +
                    gndr + hinctnta + lrscale + 
                    mbtru + seat_share +
                    f.cntry + f.inwyys, 
                  data = df.no.ch,
                  family = binomial(link = "logit"))
summary(reg2.no.ch)

reg2.no.hupl <- glm(clsprty.same ~ f.coalition.current + f.prtvtfam + f.coalition.current*f.prtvtfam +
                      agea + blgetmg + eduyrs +
                      gndr + hinctnta + lrscale + 
                      mbtru + seat_share +
                      f.cntry + f.inwyys, 
                    data = df.no.hupl,
                    family = binomial(link = "logit"))
summary(reg2.no.hupl)

reg2.no.all <- glm(clsprty.same ~ f.coalition.current + f.prtvtfam + f.coalition.current*f.prtvtfam +
                     agea + blgetmg + eduyrs +
                     gndr + hinctnta + lrscale + 
                     mbtru + seat_share +
                     f.cntry + f.inwyys, 
                   data = df.no.all,
                   family = binomial(link = "logit"))
summary(reg2.no.all)

stargazer(reg2.no.ch, reg2.no.hupl, reg2.no.all,
          keep = c("f.coalition.current",
                   "f.prtvtfam",
                   "seat_share",
                   "govt_eff",
                   "Constant"),
          dep.var.labels = c("Party Closeness"),
#          covariate.labels = c("Coalition", "Christian Dem.", "Liberal",
#                               "Social Dem.", "Communist", "Conservative", "Seat Share", "Govt Effectiveness",
#                               "Coalition $\times$ Christian Dem.", "Coalition $\times$ Liberal",
#                               "Coalition $\times$ Social Dem.", "Coalition $\times$ Communist", 
#                               "Coalition $\times$ Conservative", 
#                               "Constant"),
          notes.label = "Standard errors in parentheses",
          table.placement = "H")


